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Abstract 

We study a loop of Josephson junctions that is quenched through its critical temperature. For 
three or more junctions, symmetry breaking states can be achieved without thermal activation, in 
spite of the fact that the relaxation time is practically constant when the critical temperature is 
approached from above. The probability for these states decreases with quenching time, but the 
dependence is not allometric. For large number of junctions, cooling does not have to be fast. For 
this case, we evaluate the standard deviation of the induced flux. Our results are consistent with 
the available experimental data. 
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We consider a process in which a superconducting loop that contains n identical Joseph- 
son junctions is cooled through its critical temperature, in the absence of applied fields, and 
monitor the spontaneous generation of metastable states with persistent current. From the 
theoretical point of view, this process is enlightening, because it provides and additional 
example of a phase transition that is dominated by the time evolution of the system param- 
eters, rather than by thermal equilibrium; this is a subject that is still far from being closed, 
and is thought to be relevant both to condensed matter physics and to cosmology From 
the practical point of view, this process has significant importance, since flux trapping is a 
major obstacle for reproducible functioning of large scale ultra-high-speed superconductiv- 
ity digital applications |l| and we would like to comprehend how it depends on the system 
parameters. 



The best known theory 



or the description of dynamic phase transitions is the Kibble- 



Zurek mechanism (KZM) . It states that in these transitions the equilibrium critical 

scalings predict various aspects of the nonequilibrium dynamics of symmetry breaking, in- 
cluding the density of residual topological defects. Several numeric simulations have tested 
the predictions of the KZM, particularly the dependence of the density of defects on the 
quenching time. Some of these simulations lead to refinements of the KZM [5[ and others 
disagree with it [6fl. The KZM has also been tested in several experiments: among them we 
will mainly be interested on those performed in superconducting loops [3, |8|, |9j . 

The study of a loop of Josephson junctions is appealing, because the "rules of the game" 
are particularly simple and the system parameters can be tailored practically at will. 
This system may be neater and qualitatively different from other systems, because it can be 
trully divided into n identical subsystems, whereas in other systems the division depends on 
a continuously varying coherence length. Moreover, experimental results are already avilable 

0. 

The supercurrent through Josephson junction i is given by 



I c (t)sxnji, (1) 



where 7, is the gauge-invariant phase difference across the junction. J c vanishes above the 
critical temperature and increases when the temperature is lowered. We will consider a 
quenching process in which I c (t) grows from to I c q. This supercurrent may be interpreted 
as arising from a potential energy term — ^ i / c cos7j. This potential energy gives rise to 



metastable asymmetric states in which the system can be trapped. 

There are two typical processes for the formation of topological defects in the KZM. In one 
case defects become confined when the order parameter becomes unable to follow the change 
of the parameters of the system. Another scenario is activation due to thermal fluctuations 
close to the Ginzburg temperature. We shall see that the present system does not fit in 
either of these cases; the relaxation time does not diverge at the critical temperature and 
no activation energy is required in order to enter a metastable state. More precisely, let 
R be the resistance of each junction, L the self-inductance of the loop, let us assume that 
the capacitance is small and the resistance of the loop itself (above the critical temperature 
T c ) is much smaller than that of the junctions. Then there are two relaxation times in the 
problem: one of them is h/2eRI c , which is infinite regardless of the temperature above T c , 
and the other is L/nR, which remains constant and refers to a process that does not attempt 
to align the order parameters of neighboring segments into the same phase. In a sense, our 



problem is similar to that of decompression of He 4 from the A-line . 

Let I C Q be the maximal superconducting current through the junctions at low temperature 
and let us take R, J c o, h and 2e as units. Accordingly, the units of voltage, energy, inductance 
and time will be RI c o, 2ei?/ c o, H/2eI c o and h/2eRI C Q. 

The state of the loop will be described by the set of values {7i}. The sum of these values 
can be interpreted as minus the magnetic flux enclosed by the loop. Since we assume that 
there is no applied magnetic flux, 

n 

E^ = - LJ ' ( 2 ) 

1=1 

where I is the current around the loop. Our goal is to find the probabilities for the possible 
values of 1% after the loop has been cooled. 

The rules for the evolution of {7,} are stated in several textbooks [l0(|. The ac Josephson 
relation is 

d-fi/dt = Vi , (3) 
where V, is the voltage across junction i. The total current through junction i is 

I = Ijt + Vi + CdV/dt + I Ni , (4) 

where Iji is the supercurrent, given by Eq. (JT]), C is the capacitance of each junction, and 
I Mi is the Johnson current. We will assume that during a period of time T\ the system is 
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kept above the critical temperature and I c = 0, then, during a period of time T2 the system 
is quenched and I c grows up to J c o and finally, during a period of time r 3 , I c = I c0 . 

In most of our calculations we will assume that the capacitance is negligible. In this case, 
from Eqs. (J2J), © and (0}, 

The case C ^ will be discussed below. We integrate Eq. §5§ by Euler iterations. For 
this purpose we divide the process into short periods of time of duration At. The Johnson 
current averaged over a single period is given by 

Im = V9i/ v 7 ^ , (6) 

where gi is a random number with zero average and variance 1 and 77 = {2k B T/R) 1 / 2 , with 
ks Boltzmann's constant and T the temperature. We assume that the temperature remains 
close to the critical temperature during the entire process; accordingly, 77 will be taken as 
constant. 

We first consider the case in which the loop is cooled instantaneously, i.e., r 2 = 0. If we 
ignore the Johnson current, Eq. ([5]) is equivalent to viscosity-dominated motion of a particle 
in n-dimensional space that feels a potential energy (1/2L)(^™ =1 •y i ) 2 — I c Y17=i cos 7i- m 
this situation, {74} evolves towards a local minimum of the potential energy. During the 
first stage, I c (t) = and the only local minimum is the plane Y17=i ^ = 0- During the last 
stage, I c (t) = 1. At this stage the absolute minimum is located at the origin, 7^ = 0, but 
for sufficiently large values of L additional local minima may also exist. We ask whether 
the values {7,} could wander in the plane ^" =1 7i = and then, when the temperature is 
lowered, flow into a local minimum different from 73 = 0. For our present purpose, two 
states such that their respective values of 7^ differ by integer multiples of 2tt and ^7=1 7* ^ s 
the same for both, will be considered equivalent. 

Figure [H shows that this situation is possible for the case n = 3. The figure shows 
five evolution curves that start at the plane YH=i 1% = an d ; in the absence of thermal 
fluctuations, flow to a local minimum of the potential energy. Note that although the 
temperature is assumed to change instantaneously, Y^h=\ li builds up during a lapse of time 
of the order of h/eRI c0 . 
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FIG. 1: Several evolution lines for the case of three junctions, for L = 4 and in the absence of 
thermal fluctuations. All these curves start at the plane 71+72 + 73 = and converge to the point 
7l = 72 = -0.98, 73 = 2vr - 0.98. 

We studied the evolution including thermal fluctuations, for n = 3 and n — 4, for several 
values of rj. In all cases, we started from values of {7^} randomly located in the interval 
—71 < 7 j < 7V. In order to achieve initially an equilibrium distribution, evolution was followed 
during a period t\ with I c — 0. After that, evolution was followed during a period T3 with 
Ic — ho- in order to decide what is the "final" state, we should average over an additional 
period of time, in order to filter out thermal fluctuations. We found it easier to turn off at 
this stage the fluctuations and let the state converge to the nearest local minimum. For each 
set of values, this process was reapeated 1000 times and the probability for confinement in a 
given state was evaluated as the number of times this state was obtained, divided by 1000. 
In most cases the final state was the ground state 7 j = 0, but the first excited state was also 
reached. Due to the symmetry of the problem, these states are degenerate, i.e., the 7 j's can 
be permuted and all the signs can be inverted. For the parameters we considered, we did 
not find cases with higher excited states. 

The results are shown in Fig. [2j We avoided values of rj that might be too small to enable 
thermalization during the period 77 . The probabilities shown in the graph correspond to the 
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FIG. 2: Probability for a current-carrying metastable state as a function of the size of thermal 
fluctuations. Parameters used: n = 40000, r 2 = 0, r 3 = 10000, At = 0.1, C = 0. n = 3, L = 4; 
* n = 4, L = 4; □ n = 4, L = 2. 

total probability of reaching any of the (degenerate) excited states. As a general trend, we see 
that the probability of ending at an excited state increases with the number of junctions and 
with the normalized self-inductance. We also see that this probability is fairly independent 
of the size of thermal fluctuations, until a sufficiently large value of r\ is reached. Beyond 
this value, there is a fast decrease of this probability. 

It is reasonable to expect that the probability for the metastable state will decrease 
when the thermal energy /cbT becomes comparable to the energy barrier that confines this 
state, i.e., the difference between the energy at the saddle-point and the energy at the local 
minimum. For n = 3, L = 4, a local minimum is at 71 = 72 = —0.98,73 = 2n — 0.98, the 
corresponding saddle point is at 71 = 72 = —0.67, 73 = 7r + 0.67 and the energy difference 
is 0.25; for n = 4, L = 4, a local minimum is at 71 = 72 = 73 = —0.83,74 = 2n — 0.83, 
the saddle point is at 71 = 72 = 73 = —0.54, 74 = n + 0.54 and the energy difference is 
0.42; similarly, for n = 4, L = 2, the energy barrier is 0.13. In all cases we find that the 
probability for the metastable state decreases to about half its maximum value when the 
thermal energy is about one eighth of the barrier energy. Clearly, the precise value depends 
on r 3 ; in principle, for r 3 — > 00, the metastable state should always decay. If the thermal 
energy becomes of the order of the energy difference between the excited and the ground 
state, then the probability for the excited state will increase with temperature (equilibrium 
probability), but we are not interested in this regime. 

Let us now study the influence of the cooling time on the probability for a metastable 
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FIG. 3: Probability for a metastable state as a function of the cooling time for loops with 3 or 4 
junctions. The empty (filled) symbols correspond to I c (t) proportional to t (to t 2 ) and are fitted 
by dashed (continuous) lines. For visibility, the line for n = 3, L = 4 and I c (t) oc t 2 has been 
lowered by 0.025. Parameters used: T\ = 40000, r 2 + r 3 = 10000, At = 0.1, C = 0. Unless stated 
otherwise, rj = 0.1. 0^ = 3, L = 4; □ n = 4, L = 2; A n = 4, L = 4, rj = 0.2. 
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FIG. 4: Like Fig. O for n = 5 and n = 6. The inset shows P(t 2 ) in the range < r 2 < 10 4 . For 
visibility, the line for n = 5, L = 2 and I c (t) oc t has been raised by 0.025.0 n = 5, L = 2; A n = 5, 
L = 1; o n = 6, L = 2. 

state. At the moment that I c (t) becomes different from zero, there will be an incentive for 
leaving this state; on the other hand, as long as I c (t) is small, the confining barrier will also 
be small and the way out will be easy. We therefore expect that the trapping probability 
will decrease with r 2 . In most of our calculations we assumed that / c (^) is proportional to 
the temperature below T c and therefore increases linearly with t, but we also considered 
the case I c (t) oc t 2 , which is more realistic for strong coupling 9j. For the present purpose, 
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simulations were repeated 10 4 times. Figures [3] and H] show our results for several values of 
n, L and 77. 

The topological charge for a given state may be defined as the sum of the topological 
charges of all junctions, where the topological charge of junction i is the closest integer to 
7i/27r. For all the cases considered in Figs. EHU the topological charge was or ±1. 

The probabilities shown in Figs. [3HU are also the expectations of the absolute value of 
the topological charge. It is therefore tempting to identify this probability with the density 
of defects, and anticipate that it will decrease as a power of r 2 . However, the arguments 
that lead to the time dependence of the defect density in 4J seem to be irrelevant in the 
present case; there is no obvious way to associate the presence of topological charge to some 
primordial coherence length and, indeed, our results cannot be fitted by a power dependence. 
Denoting the probability by P, most of the curves in our results (typically for small n and 
L) can be fitted by the empiric form P(r 2 ) oc exp[— (r 2 /r )V / 2/™] in the case I c (t) oc t and 
by -P(r 2 ) oc exp [— (t 2 / t ) V 2 ^™ -1 )] in the case I c {t) oc t 2 . The characteristic time r depends 
very strongly on the number of junctions and only weakly on the size of the energy barrier 
or on the temperature. For n = 3, 16 < t < 17; for n = 4, 26 < r < 39; for n = 5, 
350 < r < 540. 

Our empirical fits suggest that for n ^> 1 the probabilities for metastable states decrease 
very slowly with the cooling time. Indeed, in the experiment that involved 214 junctions 
the distribution of permanent currents was found to be independent of the cooling time 
(up to the order of a minute). 

Part of the probabilities shown in Figs. [3H11 do not decrease at a uniform rate. Instead, 
they seem to decay in two stages. A possible explanation might be that the region in phase 
space that in the absence of thermal fluctuations would flow into a metastable state can be 
divided into two subregions, such that escape from one subregion is much easier than escape 
from the other. 

Let us now consider large values of n and L, as were encountered in the experiment. 
In this case many different final metastable states are possible, and the most significant 
experimental quantity will be the variance of the induced flux. Our results are shown in 
Fig. [51 As in the case of small values of n and L, the general trend is increase of the typical 
flux with increase of either n or L. However, these individual increases appear to saturate. 
For instance, for L < 200, the standard deviation of the flux actually decreases with n in 
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FIG. 5: Standard deviation of the induced flux, $/$o = — Si 7«> as a function of the self-inductance 
and the number of junctions. Two curves are for fixed n, and L is shown in the lower abscissa; the 
other two are for fixed L, and n is shown in the upper abscissa. The abscissas are in logarithmic 
scale. The symbols have been joined for visibility. The parameters are as in Fig. and r] = 0.1. 
Each simulation was repeated 400 times. □ n = 100; A n = 214; o L = 100; * L = 600. 

the range 100 < n < 200. 

In the experiment [2|, n = 214 and ($ 2 ) 1//2 /$o = 7.4 ± 0.7, where $ is the induced flux 
and $o the quantum of flux. Comparison with Fig. [5] indicates that 2eI c oL/h should be in 
the range between ~ 100 and ~ 150. The estimates of Ref. [3] suggest that 2eI c0 L/h ~ 600. 
Since the experimental estimate was not based on a measurement, but rather on a plausibility 
argument for the size of I c q, and the junctions are not really all identical, the agreement is 
reasonable. 

Let us finally consider the case C ^ 0. In this case we integrated Eqs. ([3]) and (HI) as a 
system of differential equations. In digital applications, a preferred value is C = h/2eI c0 R 2 , 
which provides for fast switching without oscillations. We have repeated our calculations 
for this case and for several representative values of the other parameters. We found that a 
capacity of this size has no qualitative effect. 

In summary, we have performed simulations that describe the formation of symmetry 
breaking states when a loop of Josephson junctions is quenched. Among the typical systems 
in which symmetry breaking occurs in a dynamics-dominated process, the present system 
constitutes a class of its own. Our results agree with the experiment in the case of large n 
and raise predictions for the case of small n. 
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